{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tweedie Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In insurance premium prediction problems, the total claim amount for a covered risk usually has a continuous distribution on positive values, except for the possibility of being exact zero when the claim does not occur. One standard approach in actuarial science in modeling such data is using compound Poisson models.\n",
    "\n",
    "##### Compound Poisson distribution\n",
    "\n",
    "Let $ N $ be a random variable with Poisson distribution and $ Z_1, Z_2, ... $ be independent identically distributed random variables with Gamma distribution. Define a random variable $ Z $ by\n",
    "\n",
    "$$ Z = \\begin{cases}0, & \\mbox{if}\\ N = 0\\\\Z_1 + Z_2 + ... + Z_N, & \\mbox{if}\\ N > 0\\end{cases} $$\n",
    "\n",
    "The resulting distribution of $ Z $ is called compound Poisson distribution. In the case of insurance premium prediction $ N $ referres to the number of claims, $ Z_i $ reffers to the amount of $i$-th claim. Compound Poisson distribution is a special case of Tweedie model.\n",
    "\n",
    "Log-likelihood of compound Poisson distribution can be written as\n",
    "$$ p(z) = \\frac{1}{\\phi}\\left(z \\frac{\\mu^{1-\\rho}}{1-\\rho} - \\frac{\\mu^{2-\\rho}}{2-\\rho}\\right) + a$$\n",
    "\n",
    "where $ a, \\phi, \\mu $ and $ 1 < \\rho < 2 $ are some constants.\n",
    "\n",
    "We will apply Tweedie model to an auto insurance claim dataset analyzed in Yip, Yau (2005) and Zhou, Yang, Qian (2019).\n",
    "\n",
    "##### Loading dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://cran.r-project.org/src/contrib/cplm_0.7-8.tar.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tar -xf cplm_0.7-8.tar.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install rdata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import rdata\n",
    "\n",
    "data = rdata.parser.parse_file('cplm/data/AutoClaim.RData')\n",
    "df = rdata.conversion.convert(data)['AutoClaim']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "from catboost.utils import eval_metric\n",
    "from catboost import CatBoostRegressor, Pool"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>AGE</th>\n",
       "      <th>BLUEBOOK</th>\n",
       "      <th>HOMEKIDS</th>\n",
       "      <th>KIDSDRIV</th>\n",
       "      <th>MVR_PTS</th>\n",
       "      <th>NPOLICY</th>\n",
       "      <th>RETAINED</th>\n",
       "      <th>TRAVTIME</th>\n",
       "      <th>AREA</th>\n",
       "      <th>CAR_USE</th>\n",
       "      <th>CAR_TYPE</th>\n",
       "      <th>GENDER</th>\n",
       "      <th>JOBCLASS</th>\n",
       "      <th>MAX_EDUC</th>\n",
       "      <th>MARRIED</th>\n",
       "      <th>REVOLKED</th>\n",
       "      <th>CLM_AMT5</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1019</th>\n",
       "      <td>45</td>\n",
       "      <td>14830</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>6</td>\n",
       "      <td>31</td>\n",
       "      <td>Urban</td>\n",
       "      <td>Private</td>\n",
       "      <td>Sedan</td>\n",
       "      <td>M</td>\n",
       "      <td>Professional</td>\n",
       "      <td>Masters</td>\n",
       "      <td>Yes</td>\n",
       "      <td>No</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5461</th>\n",
       "      <td>42</td>\n",
       "      <td>13770</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>14</td>\n",
       "      <td>24</td>\n",
       "      <td>Urban</td>\n",
       "      <td>Private</td>\n",
       "      <td>Sports Car</td>\n",
       "      <td>F</td>\n",
       "      <td>Professional</td>\n",
       "      <td>Bachelors</td>\n",
       "      <td>Yes</td>\n",
       "      <td>No</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7226</th>\n",
       "      <td>55</td>\n",
       "      <td>21520</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>25</td>\n",
       "      <td>Urban</td>\n",
       "      <td>Private</td>\n",
       "      <td>Van</td>\n",
       "      <td>M</td>\n",
       "      <td>Blue Collar</td>\n",
       "      <td>&lt;High School</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>6656</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6233</th>\n",
       "      <td>33</td>\n",
       "      <td>25380</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>6</td>\n",
       "      <td>27</td>\n",
       "      <td>Urban</td>\n",
       "      <td>Commercial</td>\n",
       "      <td>Panel Truck</td>\n",
       "      <td>M</td>\n",
       "      <td>Blue Collar</td>\n",
       "      <td>High School</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8215</th>\n",
       "      <td>45</td>\n",
       "      <td>22680</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>6</td>\n",
       "      <td>24</td>\n",
       "      <td>Urban</td>\n",
       "      <td>Private</td>\n",
       "      <td>Sedan</td>\n",
       "      <td>M</td>\n",
       "      <td>Professional</td>\n",
       "      <td>Masters</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>6314</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      AGE  BLUEBOOK  HOMEKIDS  KIDSDRIV  MVR_PTS  NPOLICY  RETAINED  TRAVTIME  \\\n",
       "1019   45     14830         2         0        0        3         6        31   \n",
       "5461   42     13770         3         1        0        1        14        24   \n",
       "7226   55     21520         0         0        4        1         1        25   \n",
       "6233   33     25380         0         0        0        2         6        27   \n",
       "8215   45     22680         0         0        5        1         6        24   \n",
       "\n",
       "       AREA     CAR_USE     CAR_TYPE GENDER      JOBCLASS      MAX_EDUC  \\\n",
       "1019  Urban     Private        Sedan      M  Professional       Masters   \n",
       "5461  Urban     Private   Sports Car      F  Professional     Bachelors   \n",
       "7226  Urban     Private          Van      M   Blue Collar  <High School   \n",
       "6233  Urban  Commercial  Panel Truck      M   Blue Collar   High School   \n",
       "8215  Urban     Private        Sedan      M  Professional       Masters   \n",
       "\n",
       "     MARRIED REVOLKED  CLM_AMT5  \n",
       "1019     Yes       No         0  \n",
       "5461     Yes       No         0  \n",
       "7226      No       No      6656  \n",
       "6233      No       No         0  \n",
       "8215      No       No      6314  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "features = ['AGE', 'BLUEBOOK', 'HOMEKIDS', 'KIDSDRIV', 'MVR_PTS', \n",
    "            'NPOLICY', 'RETAINED', 'TRAVTIME', 'AREA', 'CAR_USE', \n",
    "            'CAR_TYPE', 'GENDER', 'JOBCLASS', 'MAX_EDUC', 'MARRIED', \n",
    "            'REVOLKED']\n",
    "cat_features = ['AREA', 'CAR_USE', 'CAR_TYPE', 'GENDER', 'JOBCLASS', \n",
    "                'MAX_EDUC', 'MARRIED', 'REVOLKED']\n",
    "target = 'CLM_AMT5'\n",
    "df = df[features + [target]]\n",
    "df_train, df_test = train_test_split(df, random_state=0)\n",
    "\n",
    "df_train.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset contains various features describing drivers and their vehicles. The CLM_AMT5 column contains the total claim amount over last 5 years. It will be our target. The following histogram shows distribution of the target."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAWLUlEQVR4nO3df2zc9X3H8eerhB8ZpU3Cj1OURAuoUVcQg6ZWSMVUeWRNQqia/FGkVKh4LJOnNa2oFqkLq7SoUCTaidHCWjarZDNVWshoUSJGS62UU1dtCSEFEiBlMSElVjKy1iGtQaVz994f9zGck7PvbN+dffd5PaTTfb/v+3zvPm/n9Lqvv9/vOYoIzMwsD++a7gmYmVnzOPTNzDLi0Dczy4hD38wsIw59M7OMzJruCYznoosuisWLF096+zfeeIPzzz+/fhOaAdxT62jHvtxTa9i3b98vIuLiSo/N6NBfvHgxTz/99KS3LxaLdHZ21m9CM4B7ah3t2Jd7ag2Sfj7WYz68Y2aWEYe+mVlGHPpmZhlx6JuZZcShb2aWkaqhL+n9kp4tu/1K0uckzZPUJ+lQup+bxkvSvZL6Je2XtLTsubrS+EOSuhrZmJmZnalq6EfESxFxdURcDXwIeBN4FNgM7IqIJcCutA5wPbAk3bqB+wEkzQO2ANcAy4AtIx8UZmbWHBM9vLMCeDkifg6sBXpTvRdYl5bXAg9GyW5gjqT5wCqgLyIGI+Ik0AesnnIHZmZWs4l+OWs98J20XIiI4wARcVzSJam+ADhats1Aqo1VH0VSN6XfECgUChSLxQlO8R1DQ0NT2n4mck+tox37ck+tr+bQl3QO8HHgtmpDK9RinProQkQP0APQ0dERU/mm3H3bdnD3T96Y1LZH7rph0q/bSO347cF27Anasy/31PomcnjneuCnEfFaWn8tHbYh3Z9I9QFgUdl2C4Fj49TNzKxJJhL6n+SdQzsAO4GRK3C6gB1l9ZvTVTzLgVPpMNATwEpJc9MJ3JWpZmZmTVLT4R1Jvwd8FPiLsvJdwHZJG4BXgRtT/XFgDdBP6UqfWwAiYlDSHcDeNO72iBiccgdmZlazmkI/It4ELjyt9ktKV/OcPjaAjWM8z1Zg68SnaWZm9eBv5JqZZcShb2aWEYe+mVlGHPpmZhlx6JuZZcShb2aWEYe+mVlGHPpmZhlx6JuZZcShb2aWEYe+mVlGHPpmZhlx6JuZZcShb2aWEYe+mVlGHPpmZhlx6JuZZcShb2aWEYe+mVlGHPpmZhlx6JuZZaSm0Jc0R9Ijkn4m6aCkD0uaJ6lP0qF0PzeNlaR7JfVL2i9padnzdKXxhyR1NaopMzOrrNY9/a8BP4iIPwCuAg4Cm4FdEbEE2JXWAa4HlqRbN3A/gKR5wBbgGmAZsGXkg8LMzJqjauhLeg/wEeABgIj4bUS8DqwFetOwXmBdWl4LPBglu4E5kuYDq4C+iBiMiJNAH7C6rt2Ymdm4ZtUw5jLgf4B/lnQVsA+4FShExHGAiDgu6ZI0fgFwtGz7gVQbqz6KpG5KvyFQKBQoFosT6WeUwmzYdOXwpLadyus20tDQ0Iyd22S1Y0/Qnn25p9ZXS+jPApYCn42IPZK+xjuHcipRhVqMUx9diOgBegA6Ojqis7OzhilWdt+2Hdx9oJYWz3Tkpsm/biMVi0Wm8jOZidqxJ2jPvtxT66vlmP4AMBARe9L6I5Q+BF5Lh21I9yfKxi8q234hcGycupmZNUnV0I+I/waOSnp/Kq0AXgR2AiNX4HQBO9LyTuDmdBXPcuBUOgz0BLBS0tx0AndlqpmZWZPUeuzjs8A2SecAh4FbKH1gbJe0AXgVuDGNfRxYA/QDb6axRMSgpDuAvWnc7RExWJcuzMysJjWFfkQ8C3RUeGhFhbEBbBzjebYCWycyQTMzqx9/I9fMLCMOfTOzjDj0zcwy4tA3M8uIQ9/MLCMOfTOzjDj0zcwy4tA3M8uIQ9/MLCMOfTOzjDj0zcwy4tA3M8uIQ9/MLCMOfTOzjDj0zcwy4tA3M8uIQ9/MLCMOfTOzjDj0zcwy4tA3M8tITaEv6YikA5KelfR0qs2T1CfpULqfm+qSdK+kfkn7JS0te56uNP6QpK7GtGRmZmOZyJ7+H0fE1RHRkdY3A7siYgmwK60DXA8sSbdu4H4ofUgAW4BrgGXAlpEPCjMza46pHN5ZC/Sm5V5gXVn9wSjZDcyRNB9YBfRFxGBEnAT6gNVTeH0zM5ugWTWOC+CHkgL4p4joAQoRcRwgIo5LuiSNXQAcLdt2INXGqo8iqZvSbwgUCgWKxWLt3ZymMBs2XTk8qW2n8rqNNDQ0NGPnNlnt2BO0Z1/uqfXVGvrXRsSxFOx9kn42zlhVqMU49dGF0gdKD0BHR0d0dnbWOMUz3bdtB3cfqLXF0Y7cNPnXbaRischUfiYzUTv2BO3Zl3tqfTUd3omIY+n+BPAopWPyr6XDNqT7E2n4ALCobPOFwLFx6mZm1iRVQ1/S+ZIuGFkGVgLPAzuBkStwuoAdaXkncHO6imc5cCodBnoCWClpbjqBuzLVzMysSWo59lEAHpU0Mv7bEfEDSXuB7ZI2AK8CN6bxjwNrgH7gTeAWgIgYlHQHsDeNuz0iBuvWiZmZVVU19CPiMHBVhfovgRUV6gFsHOO5tgJbJz5NMzOrB38j18wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIzWHvqSzJD0j6bG0fqmkPZIOSXpY0jmpfm5a70+PLy57jttS/SVJq+rdjJmZjW8ie/q3AgfL1r8M3BMRS4CTwIZU3wCcjIj3AfekcUi6HFgPXAGsBr4h6aypTd/MzCaiptCXtBC4AfhmWhdwHfBIGtILrEvLa9M66fEVafxa4KGIeCsiXgH6gWX1aMLMzGozq8ZxXwU+D1yQ1i8EXo+I4bQ+ACxIywuAowARMSzpVBq/ANhd9pzl27xNUjfQDVAoFCgWi7X2cobCbNh05XD1gRVM5XUbaWhoaMbObbLasSdoz77cU+urGvqSPgaciIh9kjpHyhWGRpXHxtvmnUJED9AD0NHREZ2dnacPqdl923Zw94FaP9dGO3LT5F+3kYrFIlP5mcxE7dgTtGdf7qn11ZKI1wIfl7QGOA94D6U9/zmSZqW9/YXAsTR+AFgEDEiaBbwXGCyrjyjfxszMmqDqMf2IuC0iFkbEYkonYn8UETcBTwKfSMO6gB1peWdaJz3+o4iIVF+fru65FFgCPFW3TszMrKrJHfso+WvgIUlfAp4BHkj1B4BvSeqntIe/HiAiXpC0HXgRGAY2RsTvpvD6ZmY2QRMK/YgoAsW0fJgKV99ExG+AG8fY/k7gzolO0szM6sPfyDUzy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy0jV0Jd0nqSnJD0n6QVJX0z1SyXtkXRI0sOSzkn1c9N6f3p8cdlz3ZbqL0la1aimzMysslr29N8CrouIq4CrgdWSlgNfBu6JiCXASWBDGr8BOBkR7wPuSeOQdDmwHrgCWA18Q9JZ9WzGzMzGVzX0o2QorZ6dbgFcBzyS6r3AurS8Nq2THl8hSan+UES8FRGvAP3Asrp0YWZmNZlVy6C0R74PeB/wdeBl4PWIGE5DBoAFaXkBcBQgIoYlnQIuTPXdZU9bvk35a3UD3QCFQoFisTixjsoUZsOmK4erD6xgKq/bSENDQzN2bpPVjj1Be/blnlpfTaEfEb8DrpY0B3gU+EClYeleYzw2Vv301+oBegA6Ojqis7OzlilWdN+2Hdx9oKYWz3Dkpsm/biMVi0Wm8jOZidqxJ2jPvtxT65vQ1TsR8TpQBJYDcySNJOpC4FhaHgAWAaTH3wsMltcrbGNmZk1Qy9U7F6c9fCTNBv4EOAg8CXwiDesCdqTlnWmd9PiPIiJSfX26uudSYAnwVL0aMTOz6mo59jEf6E3H9d8FbI+IxyS9CDwk6UvAM8ADafwDwLck9VPaw18PEBEvSNoOvAgMAxvTYSMzM2uSqqEfEfuBD1aoH6bC1TcR8RvgxjGe607gzolP08zM6sHfyDUzy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDJSNfQlLZL0pKSDkl6QdGuqz5PUJ+lQup+b6pJ0r6R+SfslLS17rq40/pCkrsa1ZWZmldSypz8MbIqIDwDLgY2SLgc2A7siYgmwK60DXA8sSbdu4H4ofUgAW4BrgGXAlpEPCjMza46qoR8RxyPip2n518BBYAGwFuhNw3qBdWl5LfBglOwG5kiaD6wC+iJiMCJOAn3A6rp2Y2Zm45o1kcGSFgMfBPYAhYg4DqUPBkmXpGELgKNlmw2k2lj101+jm9JvCBQKBYrF4kSmOEphNmy6cnhS207ldRtpaGhoxs5tstqxJ2jPvtxT66s59CW9G/gu8LmI+JWkMYdWqMU49dGFiB6gB6CjoyM6OztrneIZ7tu2g7sPTOhz7W1Hbpr86zZSsVhkKj+Tmagde4L27Ms9tb6art6RdDalwN8WEd9L5dfSYRvS/YlUHwAWlW2+EDg2Tt3MzJqklqt3BDwAHIyIvy97aCcwcgVOF7CjrH5zuopnOXAqHQZ6AlgpaW46gbsy1czMrElqOfZxLfAp4ICkZ1Ptb4C7gO2SNgCvAjemxx4H1gD9wJvALQARMSjpDmBvGnd7RAzWpQszM6tJ1dCPiJ9Q+Xg8wIoK4wPYOMZzbQW2TmSCZmZWP/5GrplZRhz6ZmYZceibmWXEoW9mlhGHvplZRhz6ZmYZceibmWXEoW9mlhGHvplZRhz6ZmYZceibmWXEoW9mlhGHvplZRhz6ZmYZceibmWXEoW9mlhGHvplZRhz6ZmYZceibmWWklv8YPUuLN//bpLc9ctcNdZyJmVn9eE/fzCwjVUNf0lZJJyQ9X1abJ6lP0qF0PzfVJeleSf2S9ktaWrZNVxp/SFJXY9oxM7Px1LKn/y/A6tNqm4FdEbEE2JXWAa4HlqRbN3A/lD4kgC3ANcAyYMvIB4WZmTVP1dCPiB8Dg6eV1wK9abkXWFdWfzBKdgNzJM0HVgF9ETEYESeBPs78IDEzswab7IncQkQcB4iI45IuSfUFwNGycQOpNlb9DJK6Kf2WQKFQoFgsTnKKUJgNm64cnvT2kzWVOVczNDTU0OefDu3YE7RnX+6p9dX76h1VqMU49TOLET1AD0BHR0d0dnZOejL3bdvB3Qeaf4HSkZs6G/bcxWKRqfxMZqJ27Anasy/31Pomm4ivSZqf9vLnAydSfQBYVDZuIXAs1TtPqxcn+dozni/3NLOZarKXbO4ERq7A6QJ2lNVvTlfxLAdOpcNATwArJc1NJ3BXppqZmTVR1T19Sd+htJd+kaQBSlfh3AVsl7QBeBW4MQ1/HFgD9ANvArcARMSgpDuAvWnc7RFx+slhMzNrsKqhHxGfHOOhFRXGBrBxjOfZCmyd0OzMzKyu/I1cM7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLi0Dczy4hD38wsIw59M7OMOPTNzDLS/P9L0MZV7X/d2nTlMH86xhj/r1tmVo339M3MMuLQNzPLiEPfzCwjDn0zs4w49M3MMuLQNzPLSNMv2ZS0GvgacBbwzYi4q9lzMBtR7RLZqRjv8tqp8KW5NhVNDX1JZwFfBz4KDAB7Je2MiBebOY92NdUAm0qYTOW1GxWO7aqRH1TVTOXfyh9WM0Oz9/SXAf0RcRhA0kPAWsChPwNMZ5hY+5up769G7nTMxA86RUTzXkz6BLA6Iv48rX8KuCYiPlM2phvoTqvvB16awkteBPxiCtvPRO6pdbRjX+6pNfx+RFxc6YFm7+mrQm3Up05E9AA9dXkx6emI6KjHc80U7ql1tGNf7qn1NfvqnQFgUdn6QuBYk+dgZpatZof+XmCJpEslnQOsB3Y2eQ5mZtlq6uGdiBiW9BngCUqXbG6NiBca+JJ1OUw0w7in1tGOfbmnFtfUE7lmZja9/I1cM7OMOPTNzDLSlqEvabWklyT1S9o83fM5naStkk5Ier6sNk9Sn6RD6X5uqkvSvamX/ZKWlm3TlcYfktRVVv+QpANpm3slVbpUtt49LZL0pKSDkl6QdGub9HWepKckPZf6+mKqXyppT5rjw+nCBCSdm9b70+OLy57rtlR/SdKqsvq0vF8lnSXpGUmPtUNPko6k98ezkp5OtZZ+/zVERLTVjdIJ4peBy4BzgOeAy6d7XqfN8SPAUuD5stpXgM1peTPw5bS8Bvg+pe84LAf2pPo84HC6n5uW56bHngI+nLb5PnB9E3qaDyxNyxcA/wVc3gZ9CXh3Wj4b2JPmux1Yn+r/CPxlWv408I9peT3wcFq+PL0XzwUuTe/Rs6bz/Qr8FfBt4LG03tI9AUeAi06rtfT7rxG3dtzTf/tPPUTEb4GRP/UwY0TEj4HB08prgd603AusK6s/GCW7gTmS5gOrgL6IGIyIk0AfsDo99p6I+M8ovVMfLHuuhomI4xHx07T8a+AgsKAN+oqIGEqrZ6dbANcBj4zR10i/jwAr0h7hWuChiHgrIl4B+im9V6fl/SppIXAD8M20rlbvaQwt/f5rhHYM/QXA0bL1gVSb6QoRcRxKAQpckupj9TNefaBCvWnSr/8fpLRX3PJ9pcMgzwInKIXAy8DrETFcYS5vzz89fgq4kIn322hfBT4P/F9av5DW7ymAH0rap9Kfc4E2eP/VW9P/tHITVP1TDy1mrH4mWm8KSe8Gvgt8LiJ+Nc5hz5bpKyJ+B1wtaQ7wKPCBceYy0flX2vFqaF+SPgaciIh9kjpHyuPMY8b3lFwbEcckXQL0SfrZOGNb5v1Xb+24p9+qf+rhtfQrJOn+RKqP1c949YUV6g0n6WxKgb8tIr6Xyi3f14iIeB0oUjoGPEfSyE5T+Vzenn96/L2UDuVNtN9Guhb4uKQjlA69XEdpz7+VeyIijqX7E5Q+nJfRRu+/upnukwr1vlH67eUwpRNLIyeRrpjueVWY52JGn8j9O0afcPpKWr6B0Secnkr1ecArlE42zU3L89Jje9PYkRNOa5rQjygd5/zqafVW7+tiYE5ang38O/Ax4F8ZfdLz02l5I6NPem5Py1cw+qTnYUonPKf1/Qp08s6J3JbtCTgfuKBs+T+A1a3+/mvIz2q6J9CgN8AaSlePvAx8YbrnU2F+3wGOA/9LaQ9iA6VjpLuAQ+l+5I0mSv/xzMvAAaCj7Hn+jNLJs37glrJ6B/B82uYfSN+8bnBPf0Tp1939wLPptqYN+vpD4JnU1/PA36b6ZZSu5uhPYXluqp+X1vvT45eVPdcX0txfouzKj+l8vzI69Fu2pzT359LthZHXbPX3XyNu/jMMZmYZacdj+mZmNgaHvplZRhz6ZmYZceibmWXEoW9mlhGHvplZRhz6ZmYZ+X8cxUj9G+c+sQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = df[target].hist(bins=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the distribution has a point mass at 0 and is right skewed. So the use of Tweedie model is well justified.\n",
    "\n",
    "#### Tweedie loss\n",
    "\n",
    "For computational stability instead of optimizing $ \\mu $ parameter of Tweedie distribution directly, we will optimize $ \\log{\\mu} $. So the Tweedie loss is given by the following formula:\n",
    "$$L = \\sum_{i=1}^n w_i \\left(-\\frac{y_i \\exp{(F(x_i)(1-\\rho))}}{1 - \\rho} + \\frac{\\exp{(F(x_i)(2-\\rho))}}{2 - \\rho}\\right) $$\n",
    "where $ w_i $ are object weights, $y_i$ is target, $ F(x_i) $ is current object prediction, $\\rho $ is the obligatory hyperparameter variance power. Variance power must belong to the interval $ (1, 2) $. \n",
    "\n",
    "#### Fitting the model\n",
    "\n",
    "We will train two CatBoostRegressor models: one trained with Tweedie loss, the other one with RMSE loss. The features are remained unchanged, the categorical ones are specified in Pool's cat_features parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<catboost.core.CatBoostRegressor at 0x7f4e1c5d4a58>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train_pool = Pool(df_train[features], label=df_train[target],\n",
    "                  cat_features=cat_features)\n",
    "test_pool = Pool(df_test[features], label=df_test[target],\n",
    "                 cat_features=cat_features)\n",
    "\n",
    "cb_tweedie = CatBoostRegressor(loss_function='Tweedie:variance_power=1.9', n_estimators=500, silent=True)\n",
    "cb_tweedie.fit(train_pool, eval_set=test_pool)\n",
    "\n",
    "cb_rmse = CatBoostRegressor(loss_function='RMSE', n_estimators=500, silent=True)\n",
    "cb_rmse.fit(train_pool, eval_set=test_pool)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Evaluating the models\n",
    "\n",
    "We will use MSLE as evaluation metric as it works well with quantities that have exponential growth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "cb_rmse_pred = np.clip(cb_rmse.predict(test_pool), 0, None)\n",
    "cb_tweedie_pred = cb_tweedie.predict(test_pool)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSLE score:\n",
      "\ttweedie loss\t [31.676911817518143]\n",
      "\trmse loss\t [35.72356239317701]\n"
     ]
    }
   ],
   "source": [
    "print('MSLE score:')\n",
    "print('\\ttweedie loss\\t', eval_metric(df_test[target].to_numpy(), cb_tweedie_pred, 'MSLE'))\n",
    "print('\\trmse loss\\t', eval_metric(df_test[target].to_numpy(), cb_rmse_pred, 'MSLE'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the model trained with Tweedie loss outperforms the model trained with RMSE loss."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### References\n",
    "- He Zhou, Yi Yang, Wei Qian (2019), \"Tweedie Gradient Boosting for Extremely Unbalanced Zero-inflated Data\", *arxiv preprint, [arXiv:1811.10192](https://arxiv.org/abs/1811.10192)*\n",
    "- Yip, K. C. and Yau, K. K. (2005), \"On modeling claim frequency data in general insurance with extra zeros\", *Insurance: Mathematics and Economics*, 36, 153–163."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
